Meshfree method for fluctuating hydrodynamics 

Anamika Pandey a , Axel Klar a , Sudarshan Tiwari a 

a Fachbereich Mathematik, TV Kaiserslautern, D-67653 Kaiserslautern, Germany 

(N 

o 

Abstract 

In the current study a meshfree Lagrangian particle method for the Landau- 
Lifshitz Navier-Stokes (LLNS) equations is developed. The LLNS equa- 
CN tions incorporate thermal fluctuation into macroscopic hydrodynamics by 

, , the addition of white noise fluxes whose magnitudes are set by a fluctuation- 

dissipation theorem. The study focuses on capturing the correct variance 
and correlations computed at equilibrium flows, which are compared with 
c-H available theoretical values. Moreover, a numerical test for the random walk 

of standing shock wave has been considered for capturing the shock location. 
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t— i 1. Introduction 

Physical quantities which describe a macroscopic system in equilibrium 
seem to be very near to their mean value. Nevertheless, due to microscopic 

^ fluctuation, random deviation from this mean value, though small, do occur. 

t-h Thermal fluctuation are a source of noise in many system. These fluctuation 

J> play a major role in phase transitions and chemical kinetics. 

^ Investigating thermal fluctuation in the motion of fluids becomes essential 

at micro and nano scale, because of the various applications of micro and 
nano scale flow, ranging from micro-engineering to molecular biology. Micro- 
machines have a major impact on many disciplines (e.g. biology, medicine, 
optics, aerospace, and mechanical and electrical engineering) [UI2]- 
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The study of fluctuation at micro and nanoscale is particularly interesting 
when the fluid is under extreme conditions or near a hydrodynamic instabil- 
ity, e.g. the breakup of droplet in nanojet, fluid mixing in the Rayleigh- Taylor 
instability PUD]. 

The presence of thermal fluctuation becomes significant for larger Kund- 
sen number (Kn > O.Ol^The LLNS equations try to capture these thermal 
fluctuation as accurately as possible which is not possible in the case of 
Navier-Stokes equations. 

To describe the general theory of fluctuation in fluid dynamics is equiva- 
lent to setting up the "equation of motion" for fluctuating quantities. Lan- 
dau and Lifshitz introduced the appropriate additional terms in the general 
equation of fluid dynamics and gave an extended form of the Navier-Stokes 
equations. The Landau-Lifshitz Navier-Stokes equations are written as 

Ut + V.F = V.D + V.S, (1) 

where, U stands for the vector of conserved quantities, density of mass, mo- 
mentum and energy 

U-(|), (2) 

F denotes the hyperbolic flux and D denotes the diffusive flux of fluid dy- 
namic equations. F and D are given by 



(3) 



D = | r | , (4) 
r.v — q 

where v is the fluid velocity, P is the pressure and T denotes the temperature. 
t = 1] ^Vv + Vv T — -IV ■ v^j is the stress tensor, q = — kVT denotes 
the heat flux. Here r\ and k are the coefficients of viscosity and thermal 




1 Kn = j-, A denotes mean free path and L represents the characteristic length. 
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conductivity, respectively. For the given expression of r we have assumed 
the bulk viscosity to be zero. 

The expression for r and q relate these quantities to the velocity and 
temperature gradients respectively. But, in the presence of fluctuation there 
are also spontaneous local stresses and heat fluxes in the fluid, which are not 
related to velocity and temperature gradient. For these spontaneous local 
stresses tensor and heat fluxes, the LLNS equations introduce additional 
quantities in the fluid dynamic equations called stochastic flux, i.e. 



where the stochastic stress tensor (sst) S and stochastic heat flux (shf) H 
have zero mean and their covariances are given by 



where, k B is the Boltzmann's constant. 

These stochastic properties for S and H have been derived by a variety 
of approaches. Originally, these properties have been derived for equilibrium 
fluctuation [6j [TTJ [121 H3] and later the validity of the LLNS equations for 
non-equilibrium systems has been shown [14J. 

In this work a meshfree numerical scheme is developed for solving the 
LLNS equations. For simplicity, we will deal with one-dimensional system. 
The Lagrangian form of the LLNS equations for ID system in terms of prim- 
itive variables can be written as 





Cov{S t] (v,t),S kl (v',t)) = 2k BV T (s&S + SSSfr - S(r-v')5(t-t), 

(6) 

Cav(Hi(r,t), H 3 (r\t )) = 2k B KT 2 5* j 5(r - r')8(t - t), (7) 
Cov(S ij (r,t),H k (r',t')) = 0, (8) 
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u is the fluid velocity in x-direction and T is the temperature, s and h 

represent sst and shf in ID respectively. Momentum J = pu and energy 

density E = c v pT + \pu 2 is expressed in terms of p, u, T. By D/Dt we denote 

the Lagrangian derivative. We will take the above system with equation of 

state P = pRT, where R is the gas constant. We will demonstrate our result 

R 

for a mono- atomic, hard sphere gas for which R = ks/m and c v = 

7-1 

where m is molecular mass and 7(= |) is the ratio of specific heat. 

Now for the covariances of the stochastic fluxes in a ID system one obtains 
from equations ([6]), §f§ and (Jsj) 

Cov(s(x,t), s(x',t')) = ^ / dy J dy' J dz J dzCov(S xx (r,t),S xx (r,t')) 

- 8 ^ T ( 5 (x _ x ' )(5(t _ t ') (12) 



3a 



and similarly, 



Cov(h(x,t),h(x',t')) = — I dy I dy I dz I dz Cov(H x (r,t), H x (r ,t')) 



a 2 



= 2 ^^5(x-x')5(t-t). (13) 
a 

Here a represents the surface area of the system in the yz - plane. 

A number of numerical schemes have been developed for stochastic hydro- 
dynamic equations. A stochastic lattice-Boltzmann model has been devel- 
oped for simulating solid-fluid suspensions by Ladd [IS]. For modelling the 
Brownian motion of particles a similar approach has been used by Sharma 
and Patankar in [TB], where they coupled the fluctuating hydrodynamics 
equations with equations of motion for the particle. Moseler amd Landman 
[9] have used LLNS stochastic stress tensor for lubrication equations and ob- 
tain good comparison with molecular dynamics simulation for the breakup 
of nanojets. 

Serrano and Espanol [IT] have developed a thermodynamically consistent 
mesoscopic fluid particle model by casting their model into the so called 
GENERIC structure which allows to introduce thermal fluctuation. They 
describe a finite volume Lagrangian discretization of the continuum equations 
of hydrodynamics using Voronoi tessellation. A similar mesoscopic, Voronoi- 
based algorithm using the dissipative particle dynamics method has been 
used by Fabritiis et al.[18j. 
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Garcia et al.jl] have developed a simple finite difference scheme for the 
linearized LLNS equations. This scheme has been designed for specific prob- 
lems. In the context of adaptive mesh and algorithm refinement hybrid 
schemes that couple continuum and particle algorithms have been developed, 
see, for example, [21], |22j and [22] ■ Further, an algorithm refinement and 
coupling molecular dynamics simulations to numerics of stochastic hydrody- 
namic equations have been presented by Fabritiis [Mj . 

In a paper by Garcia et al. [5] CFD based schemes for stochastic PDEs 
have been developed, where numerical scheme for the full LLNS equations 
have been developed. Spatial and time correlation at equilibrium have been 
compared with theoretical values and DSMC simulations. Moreover, the 
effect of fluctuations on the shock drift has been shown and results are 
compared with DSMC simulation. The method is based on a third order, 
TVD Runge-Kutta temporal integrator (RK3) combined with a centered dis- 
cretization of hyperbolic and diffusive fluxes. This scheme also incorporates 
a specific interpolation for the required accuracy in variance. 

In this work we will present a grid free method for the LLNS equations. 
We will consider a particle method based on a least squares approach, see fHj . 
We concentrate on capturing the correct variance in equilibrium flow and 
compare the result with theoretical values. We will also show the fluctuation 
effect in the shock location for a standing shock wave, mentioned above. The 
concluding section will discuss future work. Since, the developed method is a 
grid free method and the distribution of particles (moving grid) can be quite 
arbitrary, the method is suitable for complicated geometry and multiphase 
flows, see for example, [20] for the classical case without stochastic fluxes. 

2. Numerical Method 

To extend the ideas of meshfree methods for the Navier-Stokes equa- 
tions to the LLNS equations we will consider the ID LLNS equations in 
a Lagrangian framework. To develop a meshfree framework for stochastic 
partial differential equations (SPDE) has a number of advantage for exam- 
ple, for the study of the dynamics of small particles at fluid interfaces or for 
the development of hybrid methods for coupling Monte Carlo methods for 
the Boltzmann equation and the fluctuating hydrodynamics equations. We 
mention that earlier studies have dealt with the coupling of DSMC for the 
Boltzmann equation and finite volume methods for the LLNS [25] . 
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In the Lagrangian framework an additional equation for the particle po- 
sition, i.e. 

"' (14) 



Dt 



u 



is solved together with rt9 - 11 ). Here u is the fluid velocity and x denote the 



position of particle in ID. To approximate spatial derivatives at every grid 
point is equivalent to approximate the spatial derivative at every particle 
positions. For solving the Lagrangian LLNS system given by ([9 11) together 
with (14), we fill the domain by particles^} These particles move with fluid 
velocity. Then the spatial derivative in equation ([9 11) are approximated at 
each particle position from its neighbouring particles. This reduces the given 
system of stochastic partial differential equations (SPDEs) to a system of 
stochastic ordinary differential equations (SODEs) with respect to time. 

We will use the MacCormack scheme [5] for the system of SODEs. The 
discretized form of the LLNS equations is 



Step 1. 



X: 



x? + AtuT 




+ 3* 



+ 



4 / dr] 



3 V dx 



\dx J 



du\ m fdhV 
dx) i \0x). 



dx / 



(15) 
(16) 



On \ ( ) 



dx^ 



\dx J . 
(17) 



(18) 



initially these particle are just a regular grid with equal spacing 
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Step 2. 

xT =x* + Atu* (19) 
P V= P ;-Atpt(^X (20) 



^dx 

m r 

■ * — ii* i 

p* } V dx J . ' 3 h V <9^ 2 A ' S \dx J ■ \dx J , ' V 9a; 



c„p* 1 \(9x / . 3 M \dx J J , 1 \dx 2 



Final Step. 



2 

For each of the above steps P, rj, k will be computed by 



(21) 



+(a:(s><s)>(i)> 



*r +1 = ^oc+*n (23) 

pr +1 = \ (PT + PD (24) 

= ^ (<+<*) (25) 

7; m+1 = -(T t m + T**) (26) 



P = pi?T, (27) 



5 Mk BrTi . . 

^=16?V— T - (28) 

Here, d denotes the molecular diameter, M is molecular mass, m = 0, 1, 2, . . . 
represents the time step and N is the number of particles in the domain. 

The approximation for sst and shf for each particle at any instant is 
computed as 

sT = J ^7 (^7T) K (30) 



^ m = \/l^ r(7r)2) K (31) 

where V c denotes the volume between two particle and are independent, 
identically distributed (iid), Gaussian random variables with zero mean and 
unit variance. 

The stochastic fluxes require some extra care in multi-step scheme, vari- 
ance in the stochastic flux s — (s, h) is given by, 



Var ((s m+i y 




(32) 

on neglecting the multiplicity of noise i.e. Var (s m ) = Var (s*). 

Because of the temporal averaging, the variance in the flux is reduced 
to half of its original magnitude. So, to include this observation the correct 
stochastic flux for a two step scheme will be s = \/2s instead of s [3] . 



Now we have to solve equations (15 - 29). The remaining task is to 



approximate the spatial derivatives on the right hand side of the prescribed 
equations. 

2.1. Meshfree approximation of spatial derivatives 

We will describe the least square approximation of spatial derivatives in 
ID. As mentioned earlier, in this method grid points are particle positions. 
Therefore, we have to approximate the derivatives at every particle position. 
Let f(t,x) be a scalar function at x and fi(t) its value at Xj G [0, L) for 
i = 1,2,3, iV for any instant t. Spatial derivatives of f(x) at x will 
be approximated on a set of neighbouring points. For limiting number of 
neighbouring points of x, one considers a weight function w = w(\\xi — x\\; h) 
with small compact support, where h determines the size of the support. We 
will consider a Gaussian weight function in the following form 

I II 2 \ II 

nt> . nn \ y> . 



w(x i -x;h) = { eXp P*^^J' (33) 
0, else 
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with a a positive constant, chosen as a = 6.25. h defines the neighbourhood 
radius for x. Let P(x, h) = {x{ : % — 1, 2, . . . , n} be the set of n neighbouring 
points of x in an interval of radius h. We have chosen h = 3dx, where dx is 
the initial spacing of particles. 

Suppose we want to approximate the derivatives of a function f(t, x) from 
its n neighbouring points sorted with respect to its distance from x. Consider 
Taylor's expansion of f(t,Xi) around x 



f (t, Xi) = f (t, x) + (/ (t, x)) x (Xi - x) 

1 2 
+ 2"(/ '(^ x )) x A x i- x ) + e i 



(34) 



where is the error in Taylor's expansion at the point Xi. The unknowns 
fx, fxx are computed by minimizing the error for i — 1, 2, 3, ... , n. The 
above system can be written as 



e = Ma-b 



(35) 



where, 



M = 



( 1 ( ^2 \ 

X ± - X -(X 1 -X) \ 



x 2 -x - (x 2 - xy 



\ x n X \X n X) 



J 



(36) 



a = [fx, fxxf , b = [/i - /, fi - /, . . . , /„ - ff and e = [e 1: e 2 , e 3 , . . . , e n ] T . 
For n > 2 this system will be over-determined for two unknowns f x and 

fxx- 

The unknowns a are obtained from a weighted least square method by 
minimizing the quadratic form 



J = J2 w i e i = (Ma-Pj 7 W (Ma-Pj 

i=i 



(37) 
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f w l ... \ 
w 2 ... 

\ ... w n J 
The minimization of J gives 

a = (M T iyM) -1 (M T W) b (38) 

and, finally, the required derivatives of f(t,x) as a linear combination of the 
discrete values 

3. Numerical Results 

3.1. Equilibrium 

This section gives the results of the above described method for an equi- 
librium scenario. The physical domain has been chosen such that the fluc- 
tuation in the system become significant. The parameters for the numerical 
simulation are given in Table [TJ 



Molecular diam- 


3.66 x 1(T 8 




Molecular 


mass 


6.63 x 10~ 23 


eter (Argon) 






(Argon) 






Reference mass 


1.78 x 1(T 3 




Reference 


tem- 


273 


density 






perature 






Sound speed 


30781 




Reference 


veloc- 


0.5 x 30781 








ity 






System length 


1.25 x 10" 4 




Reference 


mean 


6.26 x HT 6 


(L) 






free path 






System volume 


1.96 x 10~ 16 




Time step 




1.0 x KT 13 



Table 1: System parameter in CGS units for simulation of a dilute gas 



The initial spacing of particle will be dx = L/N, where is the total 
number of initial particles. Here, = 40 particles are considered. This 
number is not fixed during the simulation, particles are added or removed 
during the simulation. 



where 

W 
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The stability condition is found to be consistent with that suggested in 

0- 

(I u | +c s ) ^ < 1 (39) 

/4t? k \ At 1 
max I I < 40 

\3 p pc„/ Ax 2 2 

where c s is the sound speed, the bar indicates the reference value of quantities 
around which the system fluctuates. For the given reference state in Table 
[T] and the given initial spacing of particles, the time step has been chosen as 

At = i(r 13 s. 

3.1.1. Variance at Equilibrium 

The first benchmark is to capture the correct variance for fluctuations 
of the system at equilibrium. We consider a periodic domain with zero net 
flow and constant non-zero net flow. We take constant average density and 
temperature in both cases as given in Table [TJ The variance is computed from 
10 7 samples. We will calculate statistics in a global sense as given below. 

/ N a M{n) \ 

mean (p) = E(p) = ^ — ( £ £ P? ] , (41) 



n=l i=l 



Var (p) = E(p 2 ) - (E(p)) 2 

N s M(n) 



1 



En=i M(n) 



, [Pi' 2 

n=l i=l 




EEw 



AT S M(n) 

EE"."I|. («) 

n=l i=l 



Here, is the total number of samples and M(n) is the number of particle at 
time (nSkip + n)At (nSkip is the number of initial time steps for stabilizing 
the system). In the same way statistics for momentum and energy can be 
computed. 

Table [2] and [3] compare the theoretical variances which have been com- 
puted in ([7j, jl]) with measured variances from the meshfree stochastic 
scheme. 
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Variance of 

Density (p) 
Momentum (J) 
Energy (E) 


Exact value 

2.35 x lCr 8 
13.34 

2.84 x 10 10 


MacCormack 
with Meshfree 

2.11 x icr 8 

13.33 

2.68 x 10 10 


percentage error 

-10.2% 
-0.07% 
-5.6% 


Table 2: Variance in conserved quantities at equilibrium for zero net flow. 


Variance of 

Density (p) 
Momentum (J) 
Energy (E) 


Exact value 

2.35 x lCr 8 
18.91 

3.67 x 10 10 


MacCormack 
with Meshfree 

2.12 x icr 8 

19.01 

3.85 x 10 10 


percentage error 

-9.7% 

+0.05% 

+4.6% 



Table 3: Variance in conserved quantities at equilibrium for constant non-zero net flow. 



3.1.2. Time covariance at equilibrium 

We will measure the time covariance of density fluctuation. The analytical 
formula of the time covariance for the density can be written as [5] 



Cov (p(oj, t)p{ui 



t + t)) —{ ( 1 J exp {— u 2 Dtt\ H — exp {— u 2 Tt\ cos (c s ojt) 

\ 7/ 7 

3r — d 

H -uexp {— u 2 Tt\ sin(c s a;r)} * Var (p(u, t)) 

(43) 

where u = 2nn/L is the wave number, 7 is the ratio of specific heat, 

4 

D T = k/ ' pc v is the thermal diffusivity, D v = -r]/p is longitudinal kinematic 

viscosity, c s is the sound speed and T = - [D v + (7 — 1)D T ] is the sound 

attenuation coefficient. 

For numerical simulation the time covariance is estimated from the mean 
of N s samples, 

l N, 

Cov (p(u, t)p(u, t + t)) Ns = w J2 WW + r ) ( 44 ) 



N s 

s n=l 



Var(p(u,t)) = Var(R(t)) (45) 
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where, 

M(n) 

R(t) = j^j Pi M^nxjL) (46) 
We are considering the lowest wave number, i.e. n—1. 




Figure 1: Time covariance of density fluctuation for equilibrium problem on a periodic 
domain. 

In figure [TJ we compare the theoretical covariance of the density with 
the described meshfree simulation. We find a reasonable agreement of the 
results up to time (~4x 1CT 9 ), when a sound wave crossed the system, i.e. 
the comparison with the theory is only accurate for short time because of 
the finite size effect. 

3.2. N on- equilibrium 

In this last numerical test we consider a random walk of a standing shock 
wave due to spontaneous fluctuations. Our interest is the variance of the 
shock location as a function of time. The shock location is given by 

na(t) rL/2 nL/2 

a p (0 = / Phdx+ / ProIx = I p(x,t)dx (47) 

J-L/2 Ja(t) J-L/2 
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Pi ~ PR 

where, p is the instantaneous average density. 

The system parameters for the simulation of the standing shock is given 
below in the Table HI 

Mass density, temperature and velocity on both side of shocks are given 
by the Rankine-Hugoniot conditions. The boundary condition is adapted to 
the initial condition. 



System length 


5.0 x 10" 4 




Reference Mean 


6.26 x 10" 6 








free path 




System volume 


7.84 x 1(T 16 




Time step 


1.0 x 10" 13 


RHS mass den- 


1.78 x 1(T 3 




LHS mass den- 


4.07 x 10~ 3 


sity 






sity 




RHS velocity 


-61562 




LHS velocity 


-26933 


RHS sound 


30781 




LHS sound 


44373 


speed 






speed 




RHS tempera- 


273.0 




LHS tempera- 


567.0 


ture 






ture 





Table 4: System parameter in CGS units for simulation of standing shock waves, Mach 
number 2.0 

We considered two different shock strengths by considering Mach 2.0 , Mach 1.4. 
We have compared our results with the Finite Volume-third order Runge- 
Kutta (FVRK3) scheme in [5j , which already has good agreement with DSMC 
molecular simulation for prescribed Mach numbers. 
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Time (t) 



Figure 2: variance of shock location. Solid lines represent Meshfree simulation and dashed 
lines are FVRK3. 

4. Concluding Remarks 

The results of the simulation show that the Lagrangian particle scheme 
gives a good agreement with the theoretical values of variances for the con- 
served variables. This shows that the scheme is able to accurately repre- 
sent fluctuations in equilibrium flow. Further tests for the standing shock 
waves confirm that with a meshfree discretization we were able to reproduce 
stochastic drift of shock waves, as verified by comparison with the FVRK3 
scheme from [5], which has been compared with molecular simulation. 

It has already been mentioned in earlier literature that the ability of 
continuum model to accurately capture the fluctuation is very much sensitive 
to the construction of the numerical scheme. This is also true for the meshfree 
framework. Minor changes in implementation lead to significant changes in 
accuracy and behavior. In future work the above will be extended to higher 
dimension. We will further study higher dimensional incompressible flow and 
the dynamics of small particles at fluid interfaces. 
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